A microscopic instability in neutral magnetized plasmas 
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We show that in a neutral magnetized plasma there exist microscopic oscillatory modes, with 
wavelengths of the order of magnitude of the mean interparticle distance, which become unstable 
when the electron density exceeds a limit proportional to the square of the magnetic field. The model 
we consider is just a hnearization of the classical one for neutral plasmas, namely a system of elec- 
trons subjected to Coulomb interactions among themselves and with a uniform positive neutralizing 
background. This model is here dealt with as an actual many-body problem, without introducing 
any averaging over the individual particles. The expression of the density limit coincides, apart 
possibly from a numerical factor of order one, with the well-known Brillouin density limit for a 
nonneutral plasma, which has however a macroscopic origin. The density limit here found has the 
same order of magnitude as the operational density limit observed in several conventional tokamak 
devices. We finally show that, when the full electromagnetic interactions are taken into account, 
dispersion relations are obtained which for short wavelengths reduce to those obtained here in the 
purely Coulomb case, and for long wavelengths reproduce the familiar ones of MHD. 

PACS numbers: 52.35.-g, 41.20.Jb, 45.30.-l-s 



I. INTRODUCTION 

It is well known that the magnetic confinement of a 
pure electron plasma is possible only for electron densities 
Ue below the so-called Brillouin density limit ^Ij, i2j] given 
by 



,M 



n 



M 



eoB^ 



(1) 



where 77^^ = 1/2, while B, eg and rrie are the magnetic 
field, the permittivity of free space, and the electron 
mass. The condition n,, < n*^ can be equivalently ex- 
pressed 'AS rj < rj^^ in terms of the dimensionless param- 
eter 



V^ 



eoB^ 



UJZ. 



(2) 



where lUp — e^n^^je^mg, and lo^ = Bejvie are the elec- 
tron plasma frequency and cyclotron frequency respec- 
tively, — e being the electron charge. The Brillouin den- 
sity limit for a nonneutral plasma is usually derived by 
studying under a mean-field approximation the motion 
of a single generic electron of the plasma: for densities 
beyond the limit, the electrostatic repulsion due to the 
other charges cannot be counterbalanced by the Lorentz 
confining force exerted by the magnetic field. 

In the present paper we show that a density limit ex- 
actly of the form ([ij — possibly with a slightly different 
value of the critical parameter rj^ — exists for neutral 
plasmas too, inasmuch as the plasma presents an inter- 
nal dynamical instability beyond that limit. Here how- 
ever the instability can be revealed only by dealing with 
a microscopic many-body model involving the mutual 
Coulomb interaction of the individual charges, because 
it mainly concerns modes of wavelength of the order of 
the interparticle distance. 



The model we consider is just the classical one of point 
electrons obeying Newton's equations in an external mag- 
netic field, with Coulomb interactions both among them- 
selves and with a smeared-out positive neutralizing back- 
ground. Such a model is considered for example in the 
works of Bohm, Gross and Pines |3-5], and in the previ- 
ous ones of Langmuir and Tonks ^, i?,] . In the impossibil- 
ity of dealing with the general analytical solution for such 
a many-body problem, those authors introduced in the 
equations of motion some averaging with respect to the 
individual particle positions, and as a consequence were 
compelled to consider only plasma oscillations with wave- 
lengths much longer than the mean interparticle distance. 
Since, on the contrary, we are interested in the study of 
modes with short wavelengths, we have to adopt a differ- 
ent approach. We thus choose to stick with the original 
many-body problem, and to retain in the description of 
the system the microscopic coordinates of all the indi- 
vidual electrons. In order to simplify the mathematical 
equations, we then perform a linearization about an equi- 
librium configuration of the system. As already pointed 
out in a classical paper by Langmuir ^ , the equilibrium 
condition requires the electrons to lie on the sites of some 
regular lattice. It turns out that the normal mode solu- 
tions of the linearized equations can be studied analyti- 
cally, leading to dispersion relations which depend para- 
metrically on the plasma density. If one considers in par- 
ticular oscillations about a simple cubic lattice it turns 
out that, for densities larger than a limit of the form ((TJ, 
there exists a relevant fraction of modes for which the 
frequency becomes complex, and so the system becomes 
unstable. This instability concerns modes of wavelength 
of the order of the interparticle distance, and so cannot 
be revealed by the methods of Bohm, Gross and Pines, 
or by the equations of magnetohydrodynamics (MHD). 

The essential microscopic nature of the phenomenon 



discussed here also emerges from considerations of an 
energetic type. The origin of the instabihty hes in the 
fact that the equihbrium configuration here considered is 
not a minimum of the potential energy of the system. In- 
deed, it will be shown that the potential energy decreases 
for global displacements involving all the electrons of the 
system, when such displacements are described by plane 
waves with wavevectors having certain directions. In the 
absence of a magnetic field, the normal modes associated 
with such wavevectors are thus unstable. For a fixed 
density, an external magnetic field of large enough mag- 
nitude can stabilize them, but below a critical value the 
modes of shortest wavelengths along those directions be- 
come unstable. 

We will also briefly discuss the possible significance 
of this instability in connection with the density limit 
empirically encountered in the operation of tokamaks for 
fusion research on magnetic confinement, pointing out 
that in several cases the density limit found here is of the 
same order of magnitude as the experimental one. 

The equations of motion and their linearization are 
given in section |lTl together with the equation for the 
normal modes and the corresponding expression of the 
energy. In section IIIII the form of the dispersion rela- 
tions is studied in dependence of the relevant parame- 
ter r] = Loi/ojl, and the existence of the instability is 
exhibited. In section IIVI a generalization of the model 
is considered, in which the full electromagnetic interac- 
tions are introduced, including retardation and radiative 
terms. We prove that the purely Coulomb model con- 
sidered in the previous sections is recovered for short 
wavelengths, while in the long wavelength limit the dis- 
persion relations exactly coincide with those provided by 
the macroscopic equations of MHD for a low tempera- 
ture plasma. Finally, the possible physical relevance of 
the instability discussed here is briefly addressed in the 
Conclusion. 

Three appendices are devoted to the technical details 
of some calculations. They concern respectively the ef- 
fective field acting on an electron inside the plasma, the 
electrostatic energy in the equilibrium configuration of 
the plasma, and the dependence of the instability thresh- 
old on the orientation of the magnetic field. 



II. THE MODEL 

A. The equations of motion and their linearization 

Denoting by z^ the position vector of the i-ih electron, 
its equation of motion is 



TOeZj = -eej(zj, i) - ezj X B , 



(3) 



where e^ is the Coulomb field generated by all the other 
electrons and by the positive background, and B is an 
external magnetic field, which is supposed to be constant. 
In order to obtain a linearized system of equations of 
motion, we look for an equilibrium configuration for the 



electrons. If the plasma is assumed to be infinite (i.e., if 
all edge effects are neglected), it is easy to see that such 
an equilibrium configuration is given by the points of any 
arbitrary simple Bravais lattice. Naming ai, a2, a.^ the 
primitive translation vectors of the lattice, we have rig — 
\/V , where V ~ |ai •a2 xasj is the volume of the primitive 
cell, so that a — V^l'^ represents the lattice parameter. 
We denote by r„ = n\SL\+nia.i^nzSLz the position vector 
of an arbitrary point of the lattice, labelled by the triple 
of relative integers n = (rii, 712,^^3) € T? . We shaU also 
label with n the electron associated with this lattice site, 
and so we denote by Zn its position vector. Finally, we 
introduce the corresponding displacement Xn by 



At) 



-Xn(f), 



The linearized equations of motion about the chosen 
equilibrium configuration of the system can be shown, 
for ah n G Z^, to be 



TOeXn = x„ — ex„xrJ 



3eo 



47reo ^ 



■'--'n — m'^m 7 V 7 



where Dm, for me 1? ^ m 7^ 0, is a symmetric matrix 
whose elements depend on the lattice geometry, namely 



(.Umjij — 



o \Jva)i\Jra)j _ "ij 



(5) 



In order to prove Q, let us start from the general 
equation ([3]), and observe that up to first order in the 
displacements Xm of the electrons one can write 



where en and Gn are respectively the contributions to 
the field en of order zero and one in the Xm . 

It is clear that en is given by the constant electrostatic 
field generated by the background and by the electrons 
m 7^ n, when all these electrons are kept fixed at their 
equilibrium positions. Since the Bravais lattice is invari- 
ant under spatial reflections, we have Sn (rn) = 0. This 
is actually the reason why, as we said before, the points 
of the lattice are equilibrium positions for the electrons. 
Moreover, since the plasma is assumed to be globally 
neutral, the charge density of the background must be 

Pbg = "•eC. From dive^°-'(rn) = Pbg/eo = "ec/eo, as- 
suming that the lattice is isotropic it then follows that 

^(en )i/i9a;j(rn) — (^ij-nge/Seo- Hence, to first order in 
Xn we can write 



JO) 



(Zn) = ^— Xn • 

Sen 



(6) 



,(1) 



Concerning en , it is given by the sum of the Coulomb 
fields generated by all the electrons m ^ n, computed at 
order one in their displacements Xm- Such contributions 
are given by the well-known expression (see for instance 



chapter 4 of ref. ^) of the field Em generated by an 
electric dipole — eXm located at rm, that is 

Em(x) = --^f3^^y-^), (7) 

47reo V y y J 

with y = X — Tm- It follows that 

e'^\rn,t) = Y^ Em(rn) = -^:^ Yl ^n-m ' ^m , 

with D given by ([5]). From ([5]) and ([5]), equation (|3]) is 
then readily obtained. 



B. The equation for the normal modes 

In order to deal with the infinite system of linear dif- 
ferential equations @, we shall look as usual for normal 
mode solutions of the form 



x„ = C exp [i(k • r„ — ujt)] , 



(9) 



where the wavevector k, the frequency uj and the polar- 
ization vector C are constants. For such an ansatz, the 

field e„ (rn, t) becomes 

e(,i)(r„,t) = -^M(q) • Ccxp[*(k • r„ - ujt)] , (10) 

where we have introduced the real dimensionless sym- 
metric matrix 

^(q) = — V Dm exp(27riq • r^/a) , (11) 

in ^-^ 

which depends on the wavevector k and on the lattice 
parameter a only through their product 

q = ka/27r . 

Note that the series in ([TT|) converges only in an im- 
proper sense. In appendix El using techniques analogous 
to those already developed in [9| , we obtain an expression 

for M(q) given by the sum of an absolutely convergent 
series. For an isotropic lattice, this expression can be 
written as 



M.,(q) = ^-^ + 7V.,(q), 



(12) 



with 



iV,,(q) = (a/3 - 2/3)q25,, + (a - 4/3)q,g, 

+ 2(5/3 - a)qfS,, - J2 c^j(Hm, q) ■ (13) 

mT^O 

Here the sum runs over all the points of the dimensionless 
reciprocal lattice with m 7^ 0, namely 



H, 



1 



(mia2 X as -I- TO2a3 x ai -I- waai x 32) , (14) 



and the function Cy is obtained by subtracting to the 
function 



;(H,q) 



{H, + q,){Hj + qj) 
(H + q)2 



(15) 



the terms of order fc < 3 of its Taylor expansion in the 
variable q^/H about the origin. Finally, the constants a 
and /? can be numerically computed for any given geom- 
etry of the lattice using formulas (|A4|) - (IA5|) . Note that 
TrN(q) = for all q. 

The term 6ij/3 on the right-hand side of (IT^ is the 
equivalent of the so-called "Lorentz term" in the expres- 
sion of the local field inside isotropic dielectrics. Its con- 
tribution to the equation of motion exactly cancels the 
first term on the right-hand side of (|H). 

In conclusion, the normal mode ansatz (jH]) for the lin- 
earized equation of motion ^ leads to a linear equation 
for the polarization vector C, namely 



.we 



B X C-a;^C = 



-4q-C + N(q)-C 

q 



(16) 



where uip = ei/n^T/eom^ is the electron plasma frequency. 
This is an equation of the form A-C — 0, with the matrix 
A given by 



Trip 



'6r, 



g2 



^-A^.-(q) 



(17) 

where eijk denotes the completely antisymmetric tensor 
such that £123 = 1. The condition for the existence of 
nontrivial solutions is det A(q, cj) = 0. By solving this 
last equation with respect to w for a given q, one obtains 
the dispersion relation for the oscillations in our model 
of magnetized plasma. 



C. The energy 

As the normal modes discussed here have a purely elec- 
trostatic nature, it is of interest to have available an an- 
alytical expression for the electrostatic energy U of the 
system. It is easy to see that 



"^-IX- 



'n(Zn): 



(18) 



where (f>n is the electrostatic potential generated by all 
the charges of the plasma (electrons and background) ex- 
cept the electron n. The energy Uq, corresponding to the 
configuration in which all the electrons are at their equi- 
librium positions, i.e. Zj, = Vn for all n g Z'^, can be eval- 
uated by means of a suitable modification of the Ewald 
method for the calculation of the electrostatic energy of 
a ionic lattice (see for instance appendix B of [10|)- If 
N is the total number of electrons in the plasma, which 
is assumed to be so big that the surface effects can be 



neglected, we have 



Uo 



-N 



be' 



Sireoa ' 



(19) 



where the dimensionless constant b can be calculated us- 
ing the formula 



6 = 



E 



exp(-7rH^/0 



^H2 



E 

n#0 



FiV^K\/a) 
K\/a 



7 + 2v^, 



e 



(20) 



involving a positive parameter ^. In this formula the 
function F is defined as 



Fix) = - 



ds . 



(21) 



It can be proved that the right-hand side of ([^(7)) is inde- 
pendent of the value of ^, and that both sums converge 
very quickly for ^ of order unity. The calculations leading 
to (j20p are carried out in appendix [BJ 

Since e„(z„,t) — — V(?!)n(zn), it follows from (0) and 
(|S]) that up to second order in the displacements Xn we 
have 



t/ = t/o + 



6eo 



E-" 



87ren 



y^ X„Dr 



(22) 



m^n 



The total energy, which is conserved on the solutions of 
the equations of the motion, is then E — T + U, where 
T = {me/2) X^n^n i^ the kinetic energy. 

For a normal mode of the form ^ , it follows from ([^^ 
that the electrostatic energy per electron U/N is given 

by 



(23) 



Then, for the total energy per electron E/N of the normal 
mode we obtain 
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i/3-M(q) 



BTreofl 



|C|2-z^B.CxC* 
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where the last equality follows from the equation of mo- 
tion p^ . 

Equation ^F^ shows that the potential energy always 
increases when a single electron is displaced from its equi- 
librium position, all the others being kept fixed. We 
see however from (P5)l that, if for some q the matrix 



a;p(i/3 — M(q)) has a negative eigenvalue correspond- 
ing to some eigenvector C, then the potential energy is 
decreased by a simultaneous displacement of all the elec- 
trons according to the pattern described by the polariza- 
tion vector C and the wavevector q. This means that, in 
such a case, the potential energy at the equilibrium con- 
figuration does not present a minimum. This is directly 
connected with the existence of unstable modes since, ac- 
cording to the dynamical equation ((TC)) . for B = the 
squared frequency w^ of the normal modes is just given 
by an eigenvalue of the matrix a;p(l/3 — M(q)). Hence 
negative eigenvalues give rise to imaginary values of the 
frequency. We will see in the next section that negative 
eigenvalues actually exist in the case of the simple cubic 
lattice, when q is parallel to one of the three principal 
axes. 



III. UNSTABLE MODES OF OSCILLATION 



Dispersion relations and existence of unstable 
modes 



To find the dispersion relations we have to solve the lin- 
ear equation ([T5|) . This contains the function N(q) given 
by (IT3l) , which depends on the particular geometry of the 
Bravais lattice. In the present paper we limit ourselves 
to considering the easiest possible case, namely that of a 
simple cubic lattice, for which the primitive translation 
vectors are parallel to the three basic unit vectors of a 
cartesian frame: a.^ — au^, i = 1,2,3. In this case the 
reciprocal lattice is also simple cubic, and we have from 
([H)) Hm ~ m, with m e Z'^. For this lattice the geo- 
metrical constants appearing in (|13p are a = 8.9136 and 
/3 = 1.2267. Moreover, the constant b appearing in the 
expression (J19p of the electrostatic energy at equilibrium 
is 5 ^2.8373. 

Let us consider the particular case in which the 
wavevector k is parallel to a principal lattice axis, say 
U3. In such a case the matrix N(q) becomes diagonal. 
We shall denote Nii{qu3) — N22{qu3) = N{q), so that 
Nsaiqus) = ~2N{q). The graph of the function N{q) 
is displayed in Fig. [1] It is easily seen from (|16p that 
the dispersion relations, when expressed in terms of the 
quantities q and 

UJ = Uj/uJc , 

contain the single positive parameter ry — uyijui^. For 
instance, if B is also parallel to U3, then the dispersion 
relation for transversal normal modes is implicitly ex- 
pressed by the equation 



u? -uj + -qN{q) =0. 



(24) 



For these modes, the electrons move along circular orbits 
in planes orthogonal to B. 

In Fig. [2] we plot the dispersion curves obtained from 
a numerical solution of (|24p for various values of 77. This 
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FIG. 1: Graph of the function N{q). 
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FIG. 2: Dispersion curves for k and B both parallel to an axis 
of the simple cubic lattice, and circular transversal polariza- 
tion. The number next to each curve gives the corresponding 
value of the ratio rj/r)c, with r) = uj^/lu^ and -qc — 4.80. 



is the figure in which the phenomenon of microscopic 
plasma instabihties manifests itself. We see in fact that 
the curves are defined in the whole Brillouin zone Q < q < 
1/2 only for rj below a certain critical value 'qc — 4.80. 
This obviously follows from the fact that l^^ is a second 
degree equation in lo, which admits two real solutions 
only for rj < l/(4iV(q)). For rj = ijc there exists a double 
real solution for q = 1/2, hence 



1 



Vc 



4iV(l/2) 



^4.. 



If 77 > 77c, for q sufhciently near to 1/2 (i.e. for suffi- 
ciently short wavelengths) the two solutions of ([M)) are 
complex conjugate. For one of these solutions the elec- 
trons simultaneously spiral away from their equilibrium 
positions, so that the normal mode becomes unstable. 

Considering again the case in which both k and B are 
parallel to U3, from (|16l) we have for the longitudinal 



u^ = 7^{l + 2N{q)) 



(25) 



Since the right-hand side of this equality is always pos- 
itive, we see that for these modes w is real for all q and 
all 77. 

Let us now consider the case in which k is still parallel 
to U3, but B is parallel to U2. For the transversal mode 
with C parallel to B we find the relation w^ = —Lu'lN{q), 
or 



OJ 



-r]N{q) . 



(26) 



The right-hand side is in this case always negative, hence 
this equation provides two opposite imaginary values for 
CO. This means that, for all 77 > 0, there exists an unstable 
mode for which the electrons move exponentially away 
from the equilibrium positions along the direction of B. 
Furthermore, for elliptic orbits orthogonal to B we have 

S"* - [1 + (1 + N{q))Tj] u^ - N{q){l + 2N{q))rj'^ = . 

(27) 
This second degree equation in w^ always admits two real 
solutions of opposite signs, hence there exist a stable and 
an unstable mode of this type for all q and all 77 > 0. 

The dispersion curves for all the stable modes consid- 
ered above, with q parallel to a lattice axis, are reported 
in the right part of Fig. [3] for the particular value 77 = 3. 
The solid lines correspond to the cases in which B is 
parallel to q. In particular, the two lower ones refer 
to transversal modes and are derived from (|24p . They 
can thus be compared with the curves of Fig. [21 noticing 
that 7; = 3 corresponds to rj/ric = 0.625. The top solid 
line refers instead to the longitudinal modes described by 
((25)1 . Finally, the dashed line corresponds to the case in 
which B is orthogonal to q, and is derived from (P?]) . 

The left part of Fig. [3] shows the behavior of the same 
curves in the long-wavelength limit, i.e. for g ^ 1, as re- 
sulting from the full electromagnetic treatment given in 
section IIVI We will show that in this limit the disper- 
sion relations derived from our model exactly coincide 
with those provided by the equations of MHD for a zero 
temperature plasma. 



B. Estimate of the critical parameter rj 

From the cases just considered it appears that the 
threshold 77c for the onset of plasma instabilities depends 
on the angle between B and q = ka/27r. We can thus 
write 77c = 77c(cos0), with 77c(0) = 0, 77c(l) — 4.80. For 
a generic 6, since N is an increasing function of 5, the 
instability will first manifest itself at the edge q = 1/2 
of the Brillouin zone. Hence, to determine rjc (cos 0) we 
look for the solutions w of the equation det A(k, w) = 0, 
where A(k, cj) is given by formula (fT7)l for q ~ (l/2)u3. 
The critical value rjc (cos 9) is determined as the largest 77 
for which all these solutions are real (see appendix [C] for 
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FIG. 3: Dispersion curves obtained from our model for lS^/ujI — 3, at two different scales of the k axis. It is assumed that 
ujca/2-RC « 10~^, so the left graph shows the behavior for ka/2-K < 10~* of the curves of the right graph. The left graph (see 
section IIV|I is independent of the lattice structure, whereas the right graph is specific of the modes propagating along an axis 
of a simple cubic lattice. Solid lines describe modes with B parallel to k, dashed lines modes with B perpendicular to k. 




FIG. 4: Graph of the function r]c versus cos 5. 



the details of the calculation). The graph of the function 
r]c{cos9) obtained in this way is shown in Fig. S) 

Since there is no a priori correlation between the di- 
rection of B and the orientation of the cubic lattice, it 
seems reasonable to associate to our model of plasma 
the critical parameter 77*^ which is obtained by averag- 
ing the function rjc {cos 9) over the full solid angle. One 
thus finds with a numerical integration 



V 



M 



7jcicos9)d{cos9) ^ 1.40. 



Recalling the definition ^ of the parameter 77, we con- 
clude that to any given value of B one can associate a 
critical value n^ of the electronic density given by ((U , 
with 77^ = 1.40. For densities above this threshold, un- 
stable normal modes are expected to arise within the 



plasma in a significant way. 

By means of a more extensive study of the behavior 
of the matrix M(q) as a function of q, it is possible to 
show that there exist also unstable modes for which the 
wave propagates along other lattice directions. However, 
the case considered in this section, for which q is paral- 
lel to a lattice axis, is the most significant one, since for 
any such q there are two independent transversally polar- 
ized unstable modes. Moreover, the instability of these 
modes can be removed by the presence of a suitable ex- 
ternal magnetic field. This is precisely the mechanism 
which is responsible for the prediction of a density limit 
proportional to the squared magnetic field. 

Of course, we could have taken a different Bravais lat- 
tice, instead of a simple cubic one, as the equilibrium 
configuration of the electrons. In such a case, the number 
and the properties of the unstable modes would in gen- 
eral have been different, since the matrix M(q) depends 
on the specific geometry of the lattice. As a consequence, 
also the value of the critical parameter rj^ is expected 
to be dependent on the choice of the equilibrium lattice. 
However, a systematic investigation of this dependence 
falls outside the scope of the present work. 



IV. THE ELECTRODYNAMICAL EXTENSION 
OF THE MODEL 

A. The equation for the normal modes 

The characteristic feature of the present approach con- 
sists in linearizing the equations of motion of the clas- 
sical plasma model about an equilibrium configuration. 
By considering purely coulombian interactions, we have 



exhibited the existence of instabihties at short wave- 
lengths which were not revealed by other treatments of 
the model. In the present section we want to show that, 
on the other hand, the present approach exactly repro- 
duces for long wavelengths the dispersion relations which 
are usually obtained for cold plasmas by applying the 
continuum equations of MHD or the methods of refer- 
ences [2t5|. 

To this end, we note that the analytical treatment of 
plasma oscillations in the dipole approximation, which 
has been given in section |TT] for purely coulombian inter- 
actions, can be generalized in a straightforward way so as 
to make use of the full electrodynamical expression of the 
field, thus taking into account also radiative terms and 
retardation. This amounts to replacing the electrostat- 
ical expression ([7]) of the dipole field with the complete 
one (see for instance chapter 9 of ref. I^ ) 



Em(x,t) 



e 



Attcq 



„ (xm_y)y _ ^ , o (xm ■ y)y 
y y cy 



(xm • y)y x„ 



cy^ 



c?y-^ 



c^y 



(28) 



where y = x — rm. In this formula the vector Xm and 
all its time derivatives are evaluated at the retarded time 
^ret = t — y / c. The resulting equations for a normal 
mode essentially coincide with those obtained in reference 
Q , so we will here limit ourselves to briefly recalling the 
results. By summing the retarded fields Em generated 
by all the individual electrons m ^ n, we obtain in place 
of (UnD 



.(1) 



(rn,t) 



^Mq,/ -0 + - 

exp[i(k-r„-a;i)] , 



(29) 



where M(q, /) is a dimensionless symmetric matrix 
which depends on the rescaled wavevector q = ka/27r 
and the rescaled frequency / = Loa/2'iTc. A remarkable 
fact is that M(q, /) is real when its arguments q and 
/ are real. Moreover, the term (iea;^/67reoC'^)C inside 
the square brackets of ((29t is exactly cancelled by the 
term (e^/67reoc^) z n which, according to the Lorentz- 
Dirac equation (see [ll|,[l4| and chapter 17 of Q), has to 
be added to the right-hand side of the equation of mo- 
tion ([3]) in order to take into account radiation reaction. 
Comments on the profound mathematical and physical 
meaning of this cancellation are given in references [l3| 
and [S]. It has however to be noted that, in any case, the 
strength of the radiation reaction is completely negligible 
as far as the study of the dispersion relations in a plasma 
is concerned. ^^ 

It is useful to decompose the matrix M as 

M(q,/)=M'-(q,/)+N(q,/), 

where M'*(q, /) is the dominant term for long- 
wavelengths, i.e. for g ^ 1, while N(q, /) is a short- 
wavelength term, which becomes relevant only when a is 



not negligible with respect to the wavelength A — 27r/fc. 
For an isotropic lattice one obtains 



M'-(q,/)=M--(k,L.)-|--, 
where 1 is the 3x3 identity matrix, and 



(30) 



M--(k,c.) 



q2]^2 _ ^2 



P 



is the term associated with the macroscopic field E™^'^ 
inside the plasma. In the long-wavelength limit we have 
in fact 



E™^'=(x,t) 



eri£ 



eo 



.M'"^'=(k, uj) • C exp [i(k • x - cut)] 



We have finally for the short-wavelength part: 

iV,,(q, /) = [{a/3 - 2/?)q2 - {2a/3),f]S,j 

+ {a~ ^i)q^qJ + 2(5/3 - a)q^,5,. 
- X! Cjj(Hm,q,/), 

m/O 



(31) 



where the function Cy is obtained by subtracting to the 
function 



Cy (H, q, /) = 



{H,+q,){H,+q,)-fH, 
(H-l-q)2-/2 



the terms of order fc < 3 of its Taylor expansion in the 
variables (\/H and f /H about the origin. This implies 
that in the long wavelength limit (i.e. g <C 1, / ^ 1) % 
is of order four in the variables q and /. It follows that 
in this limit the leading term of N(q, /) is represented 
by the second-degree homogeneous polynomial appearing 
on the right-hand side of (1511) . 

From these results, one deduces that the linearized dy- 
namical equation for a normal mode is 



i—B X C-w^C 
me. 



= c.^[M--(k,c.)+N(q,/)].C, 



(32) 



which represents the generalization of (1161) when the full 
electrodynamic interaction is taken into account. 



B. The Coulomb limit 

It is now interesting to establish under which condi- 
tions the purely coulombian equation (J16p represents a 
good approximation of (|5^ . We first note in this re- 
spect that the proportionality factor /3 between / and 

Cj = w/wc is /3 = f/uj ~ ijJca/'2,TTC — {ucre/^Tr'^ric)^'^, 



where rj — ujp/ujf., while 



e /47reomeC is the so- 



called classical electron radius. In all experimental situ- 
ations, one always has Wc ^ c/r,, = 1.06 x 10^'^ Hz (i.e. 
B < 10^^ T), whence /3 < 1. For instance, for 77 = 3 
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and Wc w 10^^ Hz, which is a typical vahie for a tokamak, 
we find /?w5xl0~^. It follows that, for all dispersion 
curves discussed in section Hill for which uj was at most 
of order unity, one has / — l3Cb <^ 1. This implies that, 
for wavelengths not too much longer than the lattice pa- 
rameter, / is negligible with respect to q. Hence, with 
very good approximation one can operate the substitu- 
tion M(q, /) -> M(q, 0), and it is immediate to see that 
M(q, 0) is just the matrix M(q) we introduced in sec- 
tion |TT1 In particular, for / — > we have Mj'J"^'^(k, w) — >■ 

-qiqj/c^ and N(q, /) -> N(q,0) = N(q), where N(q) is 
the matrix defined by (fT3|) . It follows that the dynamical 
equation (15^ reduces to p^ in this approximation. We 
have thus verified that the formulas derived by consider- 
ing purely coulombian interactions give a fully satisfac- 
tory description of the dispersion relations discussed in 
section IIIIl 



C. The long-wavelength limit 

Let us now examine the form of the dispersion relations 
in the long- wavelength limit A 3> a. In this case g <C 1, 
so that / in general is no longer negligible with respect 
to q. On the other hand, as we have already observed, in 
this limit N(q, /) can be neglected, so that ([5^ can be 
simplified as 



i — B X C 



LO^C 



t^^M" 



= (k,a;)-C. (33) 



It is then easy to see that ([33| leads to the same disper- 
sion relations as those which are provided by the usual 
macroscopic treatment of high frequency waves in a mag- 
netized plasma. 

For instance, let us consider modes with k parallel 
to B. If C is parallel to k, then B x C = and 
M'"^'=(k,cj) • C = -C. It then follows from ^ that 
longitudinal waves (i.e. waves which involve an oscilla- 
tion of the electronic density) have frequency uj = ujp 
independently of the magnetic field B and of the wave- 
length A, for A ^ a. Hence LOp plays indeed the role 
of "plasma frequency" also in this model. Moreover, for 
transversal waves, i.e. C • k = 0, we find that the normal 
modes allowed by p3p correspond to circularly polarized 
waves with dispersion relation 



C2fc2 



1 



UJt, 



Uj(uj — Uc) 



(34) 



This result exactly coincides with that obtained in the 
approximation of MHD at zero temperature, represented 
by formula (17.35) and figures 17.4-17.5 of reference |14j . 
Similarly, let us consider the modes with k perpendic- 
ular to B and k ^ 1/a. For C parallel to B, we deduce 
from ((33)) the existence of transversal waves unaffected 
by the magnetic field, with dispersion relation 



These modes correspond to formula (16.32) and figure 
16.4 of [ij. In addition, for C • B = we have modes 
in which the electrons describe elliptical orbits in planes 
orthogonal to B. For these modes, ([55)1 provides the 
dispersion relation 



t^2(a;2_a;2) 



w^ 



a;2(a;2 



-OJi, 



^D' 



(36) 



2/ 2 



U! 



(35) 



which corresponds to formula (17.12) and figure 17.1 of 

The behavior of the dispersion curves in the long- 
wavelength limit, for B either parallel or perpendicular 
to k, is shown in the left part of Fig. |3] for oj^/ujI — 3. 
In order to compare the scales of the abscissa in the two 
graphs of this figure, recall that, for ojc ~ 10^2 Hz, we 
have ujcaj^-Kc « 5 x 10^^, so that kc/uJc — 10 corre- 
sponds to <7 = kajl-K « 5 x 10~*. This means that the 
long-wavelength region, represented in the left graph of 
Fig. [21 appears so narrow in the right graph that it be- 
comes practically invisible. It is clear however that the 
left graph displays the behavior for low k of the curves 
of the right graph. It has to be noted in particular that 
equation (|551) is satisfied for C • k = and w = inde- 
pendently of fc, provided k <C 1/a. Solutions of this type 
simply correspond to static deformations of the equilib- 
rium lattice, and represent the limit for g ^ 1 of the 
lowest solid curve in the right part of Fig. [31 



V. CONCLUSION 

The main result obtained in this paper is that the 
classical model of a neutral plasma (as constituted by 
point electrons with Coulomb interactions, moving in a 
smeared-out positive background), when linearized about 
an equilibrium position, generally presents unstable nor- 
mal modes. A relevant part of these modes is stabilized 
by an external magnetic field only for plasma densities 
below a maximal one, which is expressed by a Brillouin- 
type formula. 

A natural question is then whether our result, which es- 
sentially refers to a zero temperature situation inasmuch 
as it deals with normal modes about an equilibrium con- 
figuration, may be significant also for high temperature 
plasmas. Indeed, it is well known that disruptive insta- 
bilities occur beyond a density limit in fusion machines 
with magnetic confinement [l5|. In this connection one 
may remark that we are dealing here with an instability 
property, and that the raising of temperature tends to 
increase disorder rather than creating order. Thus the 
occurrence of an instability at zero temperature should 
imply instability at high temperatures as well. 

An indication that the instability discussed here for the 
linearized system might perhaps be of interest for fusion 
plasmas, comes from the remark that the density limit 
found here turns out in several cases to be in a fairly 
good agreement with the limit empirically encountered 



in the operation of the tokamaks for fusion research. For 
example, for a magnetic field _B = 5 T the Alcator C-Mod 
device shows a limit n*^ ~ 3.8 x 10^" m~^ [1^, whereas 
formula ([1]), with 77*^ ~ 1.4 as obtained for a simple cubic 
lattice, predicts n^^ ~ 3.4 x 10^" m^"^. Analogously, at 
B = 2 T, the DIII-D device [17| presents a density limit 
of rig^ ~ 6 X 10^^ m^-^, which has to be compared to the 
prediction n*^ ~ 5.4 x 10^^ m"'^ of formula ([1]). 

Although a B^ dependence of the density limit had 
been noticed by Granetz [l^| for the Alcator C experi- 
ment, it is generally believed that the currently available 
global set of experimental data on the density limit of 
toroidal machines is best fitted by Greenwald's empiri- 
cal scaling law [Tsl, according to which the limit is pro- 
portional to the plasma current density in the tokamak. 
However, it must be recalled that, despite the large the- 
oretical work on the subject, at the moment no widely 
accepted, first principles model for the density limit in 
tokamak devices appears to exist [l5[- Thus, the results 
here presented might provide a motivation for further ex- 
perimental investigations, in order to establish whether 
a quadratical dependence on the magnetic field may pro- 
vide a good description of the data for at least some class 
of machines. This might have relevant implications on 
the expected performances of future tokamaks. 
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Up to first order in c = Cexp(— iwt), we can write p 
p(o) _|_ p(i)^ where 



P^'H^)-^-eJ2sH^-r^) 



V 



y^ exp(iGm • x) , 



m5.^0 



p(^)(x, i) == e ^ exp(ik • r„) c • V(5^(x - r„) 

n 

= i^^c- (Gm + k) exp[i(Gm + k) • x] . 
In these formulas, the vectors 



^ 27r, , 

Gm = Tr(™ia2 X as + m2a3 x ai -f msai x a2) , 



with m G Z'^, represent the points of the reciprocal lat- 
tice. 

Using the Poisson equation A0 = — p/eo, we obtain for 
the electrostatic potential (j) the corresponding expansion 

= 0(0) +^(1), where 
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0(0) (x) 



^_Y^-Ml^2^^^ (Al) 



Veo ^^ G2 



m 



Appendix A: Calculation of the field acting on an 
electron 



We start from the expression of the charge density 
within the plasma 



p(x,t)=e^,53(x-zn(t)) 



'"(".o^^E^^pc-t^--) 



G' =G„+k 



It follows that the electric field inside the plasma is given 
in the dipole approximation by E = E*^"-' -I- E*-^', where 



E(°)(x) = -V0(°)(x) = ^ ^ ^ exp(zG^ • x) , 
^^0 ^0 ^™ 



e(i)(x,0 = -V0(i)(x,i) = -^^ ^^g;„ exp(^G:„ -x) 



(A2) 



G' =G„+k 



The contribution due to the electron n is En = En + En , where 



EW(x) 



e x-r,, 

47reo Ix - rn 



(2 



~yr~ / '^ P^ exp ip • X - rn) , 
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EW(x,i) = - exp(^k•r„) 

47reo 



x-r„ 



c-(x-r„) ■ 

■-3^ f^(x-r„) 



x-r„ 



(27r)3eo 



Using (|A2p and (|A3p . and proceeding as in [sl, we thus obtain 

eW(r„,i)= lim [E(i)(x,i) - E(,i)(x, t)! = -^ exp(ik • r„) lim /d^pe'P'^p 
x-^rn L J Keo ^^^ J P 



rf^p^^p exp[zp • (x - r„)] . 
(A3) 



Y,S\p^Ii^-q)~l 



— ^M(q) • cexp(jk • r„) 
Co 



where q = ka/27r, Hm = Gma/27r, and 

M(q) = -lim /d3pe-w=ey(p,q) 



^<53(p„H^)„l 



The above integral can be evaluated by expanding the integrand function c,;j(p, q), defined by formula (IT5|) . in powers 
of q/p about the origin. Denoting by q the term of order (q/p)'^ , the first four terms of this expansion are respectively 

„(0) _ PiPj 



.(1) 



— ftPi + qjPi ~ 2q • p^ 



P 



p. 



1 J o^ , ^^•P , P*PJ 



.(2 



.2 , 4(q-p)2 



.(3) 



1 



-y p4 



2q-pqiqj + (qiPj + qjPi) 



4(q.p) 



4^qp 



,2 , 2(q-p) 



21 



The limit for ry — > 0+ of the integral of these four terms gives a polynomial function of q whose coefficients can be 



evaluated numerically for any given lattice geometry. The remainder Cij — Cij — c. 



(0) ^(1) (2) (3) 



c^ ' of the integrand 



function is of order (q/p)^ for p — > +cxd, hence it is possible to put directly ry = before evaluating the integral. This 
procedure leads for an isotropic lattice to formulas ((TU)) and P^ - (IT51) . with 



a — lim — > 

77^0+ \ ^ 



exp (-ryH^) 277^/2 



13= lim I - y &4 cxp i-vUl.) + ^ 



(A4) 
(A5) 



Note that for an isotropic lattice one can also write 



y'cij(Hm,q)= lim V" Cy(Hm,q): 

■^^ — ' L— >+oo ^ — ^ 

m/0 m5^0,|m|<I, 



where 



r fH g^- (g^ + g»)(gj+gj-) ^ g.gj-+q'%/3-2g2j,, 2g4 



r 



Appendix B: Calculation of the electrostatic energy 
at equilibrium 



tion is 



From formula (jlSp it follows that the electrostatic en- 
ergy of our model of plasma in its equilibrium configura- 



Un = - 



Ne ,. 

lim 

2 x^o 



^(«)(x) 



Ancox 
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where (j)^ ' is the potential generated by all the charges 
of the plasma, given by (jAip . while —e/Aireox is the po- 
tential generated by the electron at the origin. 

In order to numerically compute Uq, it is convenient 
to introduce the auxiliary electrostatic potential ip gen- 
erated by an array of charge distributions, each given by 
the superposition of a point charge — e and a gaussian of 
total charge +e. This can formally be written as 



V'(x) ^ Yl ^(^ ~ ''n) ' 



where rn are the points of the Bravais lattice, and the 
function ip satisfies 



Ai/,(x) = - 



3/2 



(5^(x) - (-^ exp(-rya;' 



(Bl) 



r] being an arbitrary parameter. We have Uq — A + B, 
where 



and 



A 



B 



Ne ,. 

lim 

2 x->o 



Ne ,. 
lim 

2 X-!-0 



ip{^) + 



47reoa; 



0W(x)-^(x) 



We are now going to show that both terms A and B can 
be evaluated as the sums of rapidly convergent series. 
A standard integration of (IBll) provides 



^(x) 



Anegx 



where the function F is defined by (HD). We have 



lim 

x->0 



whence 



V^(x) + 



47reoa; 



lim 

47reo 2:^0 

2eo7r3/2 ' 



1 - F{^x) 



A = - 



Ne' 
Sneo 






J^(V^|rn|) 



(B2) 



In order to calculate B, we observe that, due to its 
lattice periodicity, ip can be Fourier expanded as 



ip{x.) = ^ -0m exp(iGin • x) , 



(B3) 



where Gm are the points of the reciprocal lattice, and 
1 



ipi, 



V 



?/j(x) exp(— iGm ■ x)(i X . 



We have 



— TT / F{^x)xdx 






e-"^ x^ dx 



AeoVT] 



and, for m 7^ 0, 
e 



i^: 



eol^lGmi Jo 
e 

e 
eo"t^|Gm|2 



F{^/rjx) sm{\Q^\x)dx 



1 



1 — exp — 



e '''^ cos(|Gm|a;)(ix 

X) 

IG^P 



4ry 



From (|XT|) and (|B3J) it then follows that 



2 



"00 + ^ I ^m - 



m#0 



eo^^lG, 



Ne' 
2^ 



--y — 



-exp 



\Gr 



4t] 



(B4) 



By putting ^ = a^ry/Tr, from (IB2p and (IB4p one finally 
obtains ^, with 6 given by pl| . 



Appendix C: Calculation ot ric{9) 

The equation det A(k, w) = 0, with A(k, w) given by 
(IT7|) and q = (l/2)u3 , can be explicitly written as 



- (1 + T])u^ - v[vZ{3Z + 2) + Z-^]i 



(CI) 



with u = ^2 = cjV^?, ^ = ^(1/2) = 0.05212, and 
^ = (1 + 3Z) cos^ 9. We see that the left-hand side of 
this equation is a third-degree polynomial in u, which we 
shall call P{u). Hence the corresponding normal modes 
will all be stable (i.e. have a real frequency) provided this 
polynomial admits three real nonnegative roots. We first 
note that P(0) < for all 77 > 0. A necessary condition 
for the existence of three positive roots is then P'(0) > 0, 
where P' denotes the derivative of P. Hence we must 
have 



•q < 



Z{3Z + 2) ' 



(C2) 



which implies in particular ^ > Z, i.e. cos9 > 
y^Z/{l + 3Z) ^ 0.2123. Whenever ^ is satisfied, it 
is easily seen that P'{u) has two positive roots, which 
we shall call ui and U2, with < ui < U2- Then P{u) 
will have three positive roots if and only if P{ui) > 
and P{u2) < 0. With some simple algebra, one sees that 
the validity of both these conditions is equivalent to the 
single inequality 

4^Z(1 + 3Z)^r]^ - [Z2(l + 3Zf - 2Z(,{18Z^ + 21Z + 5) 
+ (36^2 + 24:Z + l)C^]r]^ - 2[Z'{1 + Z) - Z(4 + 9Z)^ 
+ {l + 6Z)^2 _ 2^3j^ - Z2 -H 2Z^ - ^2 ^ . (C3) 
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It is found that, for ^ > Z, the third degree polynomial in 
r] on the left-hand side of (jCSp has a single positive root 
77, and that this root satisfies (IC2|) . Putting ric{0) — fj, 



it then follows that (|C1|) admits three real nonnegative 
roots for all r/ such that < ij < r/c (0) . 
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